function obj = locprojStateDep(varargin)
% locproj
%
%   locporj(y,x,w,h1,H,type)
%   locporj(y,x,w,h1,H,type,r,lambda)
%  
%   THETA is the object of interest. 
%   IF  LP then it contains HR coefficients on x, HR coefficients on x*state, HR coefficients on w(:,1), ..., HR coefficients on w(:,end), i.e.      HR*(2+size(w,2))
%   IF SLP then it contains K  coefficients on x, K  coefficients on x*state, HR coefficients on w(:,1), ..., HR coefficients on w(:,end), i.e. K*2+ HR*(  size(w,2))  
%   NOTE THAT THE K coefficients on x and x*state do NOT depend on the horizon. It is the B-sline that maps them into horizons, i.e. IRFs.
% 
%   BW = norminv(0.975); % 95% confidence level --> 1.96
    BW = norminv(0.95) ; % 90% confidence level --> 1.64
%   BW = norminv(0.84) ; % 68% confidence level --> 1.00 
    switch length(varargin)
        case 7
            y    = varargin{1};
            x    = varargin{2};
            x_times_s = varargin{3};
            w    = varargin{4};
            h1   = varargin{5};
            H    = varargin{6};
            type = varargin{7};

        case 9
            y    = varargin{1};
            x    = varargin{2};
            x_times_s = varargin{3};
            w    = varargin{4};
            h1   = varargin{5};
            H    = varargin{6};
            type = varargin{7};            
            r    = varargin{8}; % the difference operator of order r
            lam  = varargin{9}; % The shrinkage coefficient
        
        case 10
            y    = varargin{1};
            x    = varargin{2};
            x_times_s = varargin{3};
            w    = varargin{4} ;
            h1   = varargin{5} ;
            H    = varargin{6} ;
            type = varargin{7} ;            
            r    = varargin{8} ; % the difference operator of order r
            lam  = varargin{9} ; % The shrinkage coefficient
   undersmoothing= varargin{10}; % The confidence interval of the estimated dynamic multipliers are based on an undersmoothed estimate to reduce the extent of the bias
        otherwise
            error('wrong number of input arguments')
        
    end    

    if ~exist('undersmoothing','var')
        undersmoothing = 0.1;
    end
    obj = struct();
      
    stdLP = strcmp('reg',type); % If stdLP=1 then run Local Projections (LP) by Jorda (2005)
    T  = length(y);
    HR = H + 1 - h1;
    
    % construct the B-spline basis functions
    if ~stdLP
        B = bspline( (h1:H)' , h1 , H+1 , H+1-h1 , 3 );
        K = size( B , 2 );
    else
        K = HR;
    end
    Std_NW      =NaN(K,1);
    Std_NW_state=NaN(K,1);

    % building up the regression representation of the local projection
    idx = nan( (H+1)*T , 2 );
    Y   = nan( (H+1)*T , 1 );
    Xb  = zeros( (H+1)*T , K );
    XSb = zeros( (H+1)*T , K );
    Xc  = zeros( (H+1)*T , HR , size(w,2)+1 );
    
    w = [ ones(T,1) w ];
    
    for t = 1:T-h1
        
        idx_beg = (t-1)*HR + 1;
        idx_end = t*HR;

        idx( idx_beg:idx_end , 1 ) = t;
        idx( idx_beg:idx_end , 2 ) = h1:H;
        
        % y
        y_range = (t+h1) : min((t+H),T)';
        Y( idx_beg:idx_end ) = [ y( y_range ) ; nan(HR-length(y_range),1) ];

        % x
        if stdLP
            Xb( idx_beg:idx_end , : ) = eye(HR)*x(t);
            XSb(idx_beg:idx_end , : ) = eye(HR)*x_times_s(t);    
        else
            Xb( idx_beg:idx_end , : ) =       B*x(t)        ;   % First, we only impose smoothness on the coefficients associated with the IR, and
            XSb(idx_beg:idx_end , : ) =       B*x_times_s(t);    
        end

        % w
        for i = 1:size(w,2)
            Xc( idx_beg:idx_end , : , i ) = eye(HR)*w(t,i);   % we do not smooth the coefficients of the control variables
        end
        
    end
    
    X = [Xb XSb];
    for i = 1:size(w,2)
        X = [X Xc(:,:,i)];
    end
    
    select = isfinite(Y);  
    idx = idx(select,:);
    Y   = Y(select);
    X   = X(select,:);
    X   = sparse(X);

    % estimation
    IR       = zeros(H+1,1);
    IRstate  = zeros(H+1,1);
    if stdLP

        theta     = ( X'*X )\( X'*Y );
        IR(     (1+h1):end) = theta(  1:  K);
        IRstate((1+h1):end) = theta(K+1:2*K);
    
    else

        P = zeros( size(X,2) );
        % the matrix representation of the r-th difference operator
        D = eye(K);
        for k = 1:r 
            D = diff(D);
        end
        % The regressor matrix is rank deficient with B-splines so you need to penalize it!
        P(     1:K ,       1:K ) = D' * D;
        P((K+1:2*K),(K+1):(2*K)) = D' * D;

        theta             = ( X'*X +                lam*P )\( X'*Y );
        theta_UnderSmooth = ( X'*X + undersmoothing*lam*P )\( X'*Y );
        IR(     (1+h1):end) = B * theta(1  :  K);
        IRstate((1+h1):end) = B * theta(K+1:2*K);
    end
    
    if stdLP
    for ii=1:K
        [pos_unpack,~]=find(idx(:,2)==(ii-1+h1));
        Y_h = Y(pos_unpack);              % [     Y_h-y]
        if ii==1
            X_1 = X(pos_unpack,1:HR:end); % [full(X_h(:,1))-x]
        end
        X_h = X_1(1:size(Y_h,1),:);       % [full(X_h(:,2:end))-w(:,1:end)]
        resid_h = Y_h-X_h*theta(ii:HR:end);
        nwerr = ConstructStdErrLP(X_h, resid_h, 2*HR+1);
        Std_NW(ii)       = nwerr(1);
        Std_NW_state(ii) = nwerr(2);
    end
    else
        resid_h = Y-X*theta_UnderSmooth; % the confidence interval of the estimated dynamic multipliers are based on an undersmoothed estimate of THETA to reduce the extent of the bias.
        nwerr = ConstructStdErrSLP(X, resid_h, 2*H+1,lam*undersmoothing,P);
        Std_NW       = nwerr(1  :  K);
        Std_NW_state = nwerr(K+1:2*K);
    end
    
    % pack everything up
    obj.T   = T;
    obj.HR  = HR;
    obj.K   = K;
    
    if ~stdLP
        obj.B   = B;
    end
    
    obj.idx    = idx;
    obj.Y      = Y;
    obj.X      = X;
    obj.theta  = theta;
    obj.IR     = IR;
    obj.IRstate= IRstate;
    if stdLP
        obj.Interval_up      =IR      + [NaN(h1,1);    ( Std_NW      )*BW];        obj.Interval_down      =IR      + [NaN(h1,1);    (-Std_NW      )*BW];
        obj.Interval_up_state=IRstate + [NaN(h1,1);    ( Std_NW_state)*BW];        obj.Interval_down_state=IRstate + [NaN(h1,1);    (-Std_NW_state)*BW];
    else
        obj.Interval_up      =IR      + [NaN(h1,1);  B*( Std_NW      )*BW];        obj.Interval_down      =IR      + [NaN(h1,1);  B*(-Std_NW      )*BW];
        obj.Interval_up_state=IRstate + [NaN(h1,1);  B*( Std_NW_state)*BW];        obj.Interval_down_state=IRstate + [NaN(h1,1);  B*(-Std_NW_state)*BW];
    end

    % debug stuff
    % Bs is for display / debugging pourposes only
    % Bs = bspline( (h1:0.1:H)' , h1 , H+1 , H+1-h1 , 3 );
    % obj.Bs = Bs;
    % plot((h1:0.1:H)',Bs); hold on; 
    % plot((h1:0.1:H)',Bs(:,9)) since the index nodes are: -2,-1,0,1,2,3,4,5,6 (which is in position 9)
end

function B = bspline(x, xl, xr, ndx, bdeg)
    dx = (xr - xl) / ndx;
    t = xl + dx * [-bdeg:ndx-1];
    T = (0 * x + 1) * t;
    X = x * (0 * t + 1);
    P = (X - T) / dx;
    B = (T <= X) & (X < (T + dx));
    r = [2:length(t) 1];
    for k = 1:bdeg
        B = (P .* B + (k + 1 - P) .* B(:, r)) / k;
    end
end

%-------------------------------------------------------------------------%
function nwerr = ConstructStdErrLP(x, resid, nlag)
[nobs, nvar] = size(x);

xpxi = (x'*x)\eye(size(x,2));

% perform Newey-West correction
emat = [];
for i=1:nvar
    emat = [emat
        resid'];
end

hhat=emat.*x';
G=zeros(nvar,nvar); w=zeros(2*nlag+1,1);
a=0;

while a~=nlag+1
    ga=zeros(nvar,nvar);
    w(nlag+1+a,1)=(nlag+1-a)/(nlag+1);
    za=hhat(:,(a+1):nobs)*hhat(:,1:nobs-a)';
    if a==0
        ga=ga+za;
    else
        ga=ga+za+za';
    end
    G=G+w(nlag+1+a,1)*ga;
    a=a+1;
end % end of while

V=xpxi*G*xpxi;
nwerr= sqrt(diag(V));
end

%-------------------------------------------------------------------------%
function nwerr = ConstructStdErrSLP(x, resid, nlag, lambda, P)
[nobs, nvar] = size(x);

% xpxi = (x'*x)\eye(size(x,2));
  xpxi = (x'*x + lambda*P)\eye(size(x,2));

% perform Newey-West correction
emat = [];
for i=1:nvar
    emat = [emat
        resid'];
end

hhat=emat.*x';
G=zeros(nvar,nvar); w=zeros(2*nlag+1,1);
a=0;

while a~=nlag+1
    ga=zeros(nvar,nvar);
    w(nlag+1+a,1)=(nlag+1-a)/(nlag+1);
    za=hhat(:,(a+1):nobs)*hhat(:,1:nobs-a)';
    if a==0
        ga=ga+za;
    else
        ga=ga+za+za';
    end
    G=G+w(nlag+1+a,1)*ga;
    a=a+1;
end % end of while

V=xpxi*G*xpxi;
nwerr= sqrt(diag(V));
end